{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic optimisation example\n",
    "\n",
    "This example optimises a single node with only a wind, solar, gas and lignite generator under uncertainty about the gas price.\n",
    "\n",
    "In Stage 1 decisions are made about capacities of the generators while the gas price is unknown.\n",
    "\n",
    "\n",
    "First we solve assuming knowledge about the gas price, then stochastically according to the probability distribution of gas prices.\n",
    "\n",
    "We then show that the average total cost of the system of the stochastically optimised capacities is lower than the means of the solutions from the deterministically determined capacities.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Required data\n",
    "\n",
    "For this example, we need solar and wind generation time-series. For convenience, we will be fetching the time-series data directly from the renewables.ninja server. An arbitrary example of Germany's data is retrieved. \n",
    "\n",
    "The fetched files: \n",
    "- PV (1985-2016, SARAH) (6.37 MB)\n",
    "- Wind (Current fleet, onshore/offshore separate, MERRA-2) (13.93 MB)\n",
    "\n",
    "See: https://www.renewables.ninja/ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "from io import StringIO\n",
    "\n",
    "import pypsa, pandas as pd, numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from xarray import DataArray\n",
    "from pypsa.optimization.common import reindex\n",
    "from pypsa.descriptors import get_switchable_as_dense as get_as_dense\n",
    "from pypsa.descriptors import nominal_attrs\n",
    "from pypsa.descriptors import (\n",
    "    additional_linkports,\n",
    "    expand_series,\n",
    "    get_activity_mask,\n",
    "    get_bounds_pu,\n",
    ")\n",
    "from linopy.expressions import LinearExpression, merge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Retrieve PV & Wind data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "urls = {\n",
    "    \"solar_pu\": \"https://www.renewables.ninja/country_downloads/DE/ninja_pv_country_DE_sarah_corrected.csv\",\n",
    "    \"wind_pu\": \"https://www.renewables.ninja/country_downloads/DE/ninja_wind_country_DE_current-merra-2_corrected.csv\",\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fetch_timeseries_data(url):\n",
    "    \"\"\"Fetch the timeseries data from the renewable.ninja server\"\"\"\n",
    "\n",
    "    response = requests.get(url)\n",
    "    response.raise_for_status()  # Raise an error for bad responses\n",
    "\n",
    "    return pd.read_csv(\n",
    "        StringIO(response.text), skiprows=2, parse_dates=[\"time\"], index_col=\"time\"\n",
    "    )[\"national\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar_pu = fetch_timeseries_data(urls[\"solar_pu\"])\n",
    "wind_pu = fetch_timeseries_data(urls[\"wind_pu\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Major settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenarios = [\"low\", \"med\", \"high\"]\n",
    "\n",
    "# this just determines the default scenario when building stochastic model\n",
    "base_scenario = \"low\"\n",
    "\n",
    "# in EUR/MWh_th\n",
    "gas_prices = {\"low\": 40, \"med\": 70, \"high\": 100}\n",
    "\n",
    "probability = {\"low\": 0.4, \"med\": 0.3, \"high\": 0.3}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# years for weather data (solar is 1985-2015 inclusive, wind is 1980-2019)\n",
    "year_start = 2015\n",
    "year_end = 2015\n",
    "\n",
    "# 1 is hourly, 3 is 3-hourly\n",
    "frequency = 3\n",
    "\n",
    "# Fixed load in MW\n",
    "load = 1\n",
    "\n",
    "# https://github.com/ERGO-Code/HiGHS\n",
    "solver_name = \"highs\"\n",
    "\n",
    "cts = [\"DE\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assumptions = pd.DataFrame(\n",
    "    columns=[\"FOM\", \"discount rate\", \"efficiency\", \"investment\", \"lifetime\"],\n",
    "    index=[\"default\", \"onshore wind\", \"utility solar PV\", \"gas CCGT\", \"lignite\"],\n",
    ")\n",
    "\n",
    "assumptions.at[\"default\", \"FOM\"] = 3.0\n",
    "assumptions.at[\"default\", \"discount rate\"] = 0.03\n",
    "assumptions.at[\"default\", \"lifetime\"] = 25\n",
    "\n",
    "assumptions.at[\"onshore wind\", \"investment\"] = 2e6\n",
    "assumptions.at[\"utility solar PV\", \"investment\"] = 10e5\n",
    "assumptions.at[\"gas CCGT\", \"investment\"] = 7e5\n",
    "assumptions.at[\"gas CCGT\", \"efficiency\"] = 0.6\n",
    "\n",
    "assumptions.at[\"lignite\", \"investment\"] = 15e5\n",
    "assumptions.at[\"lignite\", \"efficiency\"] = 0.3\n",
    "\n",
    "# fill defaults\n",
    "assumptions = assumptions.fillna(\n",
    "    {\n",
    "        \"FOM\": assumptions.at[\"default\", \"FOM\"],\n",
    "        \"discount rate\": assumptions.at[\"default\", \"discount rate\"],\n",
    "        \"lifetime\": assumptions.at[\"default\", \"lifetime\"],\n",
    "    }\n",
    ")\n",
    "\n",
    "\n",
    "def annuity(lifetime, rate):\n",
    "    if rate == 0.0:\n",
    "        return 1 / lifetime\n",
    "    else:\n",
    "        return rate / (1.0 - 1.0 / (1.0 + rate) ** lifetime)\n",
    "\n",
    "\n",
    "# annualise investment costs, add FOM\n",
    "assumptions[\"fixed\"] = [\n",
    "    (annuity(v[\"lifetime\"], v[\"discount rate\"]) + v[\"FOM\"] / 100.0) * v[\"investment\"]\n",
    "    for i, v in assumptions.iterrows()\n",
    "]\n",
    "\n",
    "assumptions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Required functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare base network (without stochastic optimisation)\n",
    "def prepare_network(cts, gas_price):\n",
    "    network = pypsa.Network()\n",
    "\n",
    "    snapshots = pd.date_range(\n",
    "        \"{}-01-01\".format(year_start),\n",
    "        \"{}-12-31 23:00\".format(year_end),\n",
    "        freq=str(frequency) + \"H\",\n",
    "    )\n",
    "\n",
    "    network.set_snapshots(snapshots)\n",
    "\n",
    "    network.snapshot_weightings = pd.Series(float(frequency), index=network.snapshots)\n",
    "\n",
    "    for ct in cts:\n",
    "        network.add(\"Bus\", ct)\n",
    "        network.add(\"Load\", ct, bus=ct, p_set=load)\n",
    "\n",
    "        network.add(\n",
    "            \"Generator\",\n",
    "            ct + \" solar\",\n",
    "            bus=ct,\n",
    "            p_max_pu=solar_pu,\n",
    "            p_nom_extendable=True,\n",
    "            marginal_cost=0.01,  # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind\n",
    "            capital_cost=assumptions.at[\"utility solar PV\", \"fixed\"],\n",
    "        )\n",
    "\n",
    "        network.add(\n",
    "            \"Generator\",\n",
    "            ct + \" wind\",\n",
    "            bus=ct,\n",
    "            p_max_pu=wind_pu,\n",
    "            p_nom_extendable=True,\n",
    "            marginal_cost=0.02,  # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind\n",
    "            capital_cost=assumptions.at[\"onshore wind\", \"fixed\"],\n",
    "        )\n",
    "\n",
    "        network.add(\n",
    "            \"Generator\",\n",
    "            ct + \" gas\",\n",
    "            bus=ct,\n",
    "            p_nom_extendable=True,\n",
    "            efficiency=assumptions.at[\"gas CCGT\", \"efficiency\"],\n",
    "            marginal_cost=gas_price / assumptions.at[\"gas CCGT\", \"efficiency\"],\n",
    "            capital_cost=assumptions.at[\"gas CCGT\", \"fixed\"],\n",
    "        )\n",
    "\n",
    "        network.add(\n",
    "            \"Generator\",\n",
    "            ct + \" lignite\",\n",
    "            bus=ct,\n",
    "            p_nom_extendable=True,\n",
    "            efficiency=assumptions.at[\"lignite\", \"efficiency\"],\n",
    "            marginal_cost=150,\n",
    "            capital_cost=assumptions.at[\"gas CCGT\", \"fixed\"],\n",
    "        )\n",
    "\n",
    "    return network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add additional operational scenarios to the base model\n",
    "def prepare_stochastic_model(n):\n",
    "    m = n.optimize.create_model()\n",
    "\n",
    "    nonbase_scenarios = scenarios.copy()\n",
    "    nonbase_scenarios.remove(base_scenario)\n",
    "\n",
    "    # we only have generators in this example, which simplifies things\n",
    "    c = \"Generator\"\n",
    "    sns = n.snapshots\n",
    "    attr = \"p\"\n",
    "    active = None\n",
    "    column = \"bus\"\n",
    "    sign = 1\n",
    "    ext_i = n.get_extendable_i(c)\n",
    "    min_pu, max_pu = map(DataArray, get_bounds_pu(n, c, sns, ext_i, attr))\n",
    "    capacity = n.model[f\"{c}-{nominal_attrs[c]}\"]\n",
    "\n",
    "    for scenario in nonbase_scenarios:\n",
    "        # add extra operational variables for each non-base scenario\n",
    "        dispatch = m.add_variables(\n",
    "            coords=m[\"Generator-p\"].coords, name=f\"Generator-p-{scenario}\"\n",
    "        )\n",
    "        dispatch = reindex(dispatch, c, ext_i)\n",
    "\n",
    "        # add dispatch constraints\n",
    "        lhs = dispatch - max_pu * capacity  # instead of the tuple formulation\n",
    "        m.add_constraints(lhs, \"<=\", 0, f\"{c}-ext-{attr}-upper-{scenario}\", active)\n",
    "\n",
    "        lhs = dispatch - min_pu * capacity\n",
    "        m.add_constraints(lhs, \">=\", 0, f\"{c}-ext-{attr}-lower-{scenario}\", active)\n",
    "\n",
    "        # add nodal balance constraints\n",
    "        exprs = []\n",
    "        expr = DataArray(sign) * m[f\"{c}-{attr}-{scenario}\"]\n",
    "        buses = n.df(c)[column].rename(\"Bus\")\n",
    "        expr = expr.groupby(\n",
    "            buses.to_xarray()\n",
    "        ).sum()  # for linopy >=0.2, see breaking changes log\n",
    "        exprs.append(expr)\n",
    "        lhs = merge(exprs).reindex(\n",
    "            Bus=n.buses.index, fill_value=LinearExpression.fill_value\n",
    "        )\n",
    "        rhs = (\n",
    "            (-get_as_dense(n, \"Load\", \"p_set\", sns) * n.loads.sign)\n",
    "            .groupby(n.loads.bus, axis=1)\n",
    "            .sum()\n",
    "            .reindex(columns=n.buses.index, fill_value=0)\n",
    "        )\n",
    "        rhs.index.name = \"snapshot\"\n",
    "        rhs = DataArray(rhs)\n",
    "        mask = None\n",
    "        m.add_constraints(lhs, \"=\", rhs, f\"Bus-nodal_balance-{scenario}\", mask=mask)\n",
    "\n",
    "    # define the new objective\n",
    "\n",
    "    objective = []\n",
    "    weighting = n.snapshot_weightings.objective\n",
    "    weighting = weighting.loc[sns]\n",
    "    cost = (\n",
    "        get_as_dense(n, c, \"marginal_cost\", sns)\n",
    "        .loc[:, lambda ds: (ds != 0).all()]\n",
    "        .mul(weighting, axis=0)\n",
    "    )\n",
    "\n",
    "    for scenario in scenarios:\n",
    "        cost_modified = cost.copy()\n",
    "\n",
    "        if scenario == base_scenario:\n",
    "            name = f\"{c}-{attr}\"\n",
    "        else:\n",
    "            name = f\"{c}-{attr}-{scenario}\"\n",
    "            cost_modified[\"DE gas\"] = (\n",
    "                cost_modified[\"DE gas\"]\n",
    "                * gas_prices[scenario]\n",
    "                / gas_prices[base_scenario]\n",
    "            )\n",
    "\n",
    "        operation = m[name].sel({\"snapshot\": sns, c: cost.columns})\n",
    "        objective.append((operation * (probability[scenario] * cost_modified)).sum())\n",
    "\n",
    "    ext_i = n.get_extendable_i(c)\n",
    "    cost = n.df(c)[\"capital_cost\"][ext_i]\n",
    "    objective.append((capacity * cost).sum())\n",
    "\n",
    "    m.objective = merge(objective)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check that network is created correctly:\n",
    "# gas_price = 30\n",
    "# n = prepare_network(cts,gas_price)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First solve capacities for each scenario deterministically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = None\n",
    "\n",
    "for scenario in scenarios:\n",
    "    gas_price = gas_prices[scenario]\n",
    "\n",
    "    n = prepare_network(cts, gas_price)\n",
    "\n",
    "    n.optimize(solver_name=solver_name)\n",
    "\n",
    "    if results is None:\n",
    "        results = pd.DataFrame(columns=n.generators.index)\n",
    "        results.index.name = \"scenario\"\n",
    "\n",
    "    results.loc[scenario] = n.generators.p_nom_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now solve the full problem stochastically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gas_price = gas_prices[base_scenario]\n",
    "\n",
    "n = prepare_network(cts, gas_price)\n",
    "\n",
    "prepare_stochastic_model(n)\n",
    "\n",
    "n.optimize.solve_model(solver_name=solver_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.loc[\"stochastic\"] = n.generators.p_nom_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now test each set of capacities against realisations of the gas price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for scenario in scenarios:\n",
    "    gas_price = gas_prices[scenario]\n",
    "    n = prepare_network(cts, gas_price)\n",
    "    n.generators.p_nom_extendable = False\n",
    "\n",
    "    for capacity_scenario in results.index:\n",
    "        n.generators.p_nom = results.loc[capacity_scenario, n.generators.index]\n",
    "\n",
    "        print(n.generators.p_nom)\n",
    "\n",
    "        n.optimize(solver_name=solver_name)\n",
    "\n",
    "        results.at[capacity_scenario, f\"gas-p-{scenario}\"] = n.generators_t.p[\n",
    "            \"DE gas\"\n",
    "        ].sum()\n",
    "        results.at[capacity_scenario, f\"lignite-p-{scenario}\"] = n.generators_t.p[\n",
    "            \"DE lignite\"\n",
    "        ].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for capacity_scenario in results.index:\n",
    "    for g in n.generators.index:\n",
    "        results.at[capacity_scenario, f\"{g} CC\"] = (\n",
    "            results.at[capacity_scenario, g] * n.generators.at[g, \"capital_cost\"]\n",
    "        )\n",
    "\n",
    "    for scenario in scenarios:\n",
    "        results.at[capacity_scenario, f\"DE gas-{scenario} MC\"] = (\n",
    "            n.snapshot_weightings.objective.mean()\n",
    "            * gas_prices[scenario]\n",
    "            / n.generators.at[\"DE gas\", \"efficiency\"]\n",
    "            * results.at[capacity_scenario, f\"gas-p-{scenario}\"]\n",
    "        )\n",
    "        results.at[capacity_scenario, f\"DE lignite-{scenario} MC\"] = (\n",
    "            n.snapshot_weightings.objective.mean()\n",
    "            * n.generators.at[\"DE lignite\", \"marginal_cost\"]\n",
    "            * results.at[capacity_scenario, f\"lignite-p-{scenario}\"]\n",
    "        )\n",
    "\n",
    "    results.at[capacity_scenario, \"DE gas-mean MC\"] = sum(\n",
    "        [\n",
    "            probability[scenario]\n",
    "            * results.at[capacity_scenario, f\"DE gas-{scenario} MC\"]\n",
    "            for scenario in scenarios\n",
    "        ]\n",
    "    )\n",
    "    results.at[capacity_scenario, \"DE lignite-mean MC\"] = sum(\n",
    "        [\n",
    "            probability[scenario]\n",
    "            * results.at[capacity_scenario, f\"DE lignite-{scenario} MC\"]\n",
    "            for scenario in scenarios\n",
    "        ]\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, len(results.index), figsize=(len(results.index) * 4, 4))\n",
    "\n",
    "colors = {\n",
    "    \"wind\": \"b\",\n",
    "    \"solar\": \"y\",\n",
    "    \"lignite\": \"black\",\n",
    "    \"gas\": \"brown\",\n",
    "    \"gas MC\": \"orange\",\n",
    "    \"lignite MC\": \"gray\",\n",
    "}\n",
    "\n",
    "# fig.suptitle('Horizontally stacked subplots')\n",
    "\n",
    "for i, capacity_scenario in enumerate(results.index):\n",
    "    ax = axes[i]\n",
    "\n",
    "    df = pd.DataFrame(index=scenarios + [\"mean\"])\n",
    "\n",
    "    for tech in [\"solar\", \"wind\", \"gas\", \"lignite\"]:\n",
    "        df[tech] = results.at[capacity_scenario, f\"DE {tech} CC\"]\n",
    "\n",
    "    for scenario in scenarios + [\"mean\"]:\n",
    "        df.at[scenario, \"gas MC\"] = results.at[\n",
    "            capacity_scenario, f\"DE gas-{scenario} MC\"\n",
    "        ]\n",
    "        df.at[scenario, \"lignite MC\"] = results.at[\n",
    "            capacity_scenario, f\"DE lignite-{scenario} MC\"\n",
    "        ]\n",
    "\n",
    "    df.plot(kind=\"bar\", stacked=True, ax=ax, color=colors)\n",
    "\n",
    "    ax.set_title(f\"capacity scenario {capacity_scenario}\")\n",
    "\n",
    "    ax.legend(loc=\"upper left\")\n",
    "\n",
    "    ax.set_ylim([0, 2.5e6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n",
    "\n",
    "df = (\n",
    "    results[\n",
    "        [\n",
    "            \"DE solar CC\",\n",
    "            \"DE wind CC\",\n",
    "            \"DE gas CC\",\n",
    "            \"DE lignite CC\",\n",
    "            \"DE gas-mean MC\",\n",
    "            \"DE lignite-mean MC\",\n",
    "        ]\n",
    "    ]\n",
    "    .rename(columns=lambda x: x[3:-3])\n",
    "    .rename(columns={\"gas-mean\": \"gas MC\", \"lignite-mean\": \"lignite MC\"})\n",
    ")\n",
    "\n",
    "df.plot(kind=\"bar\", stacked=True, ax=ax, color=colors)\n",
    "\n",
    "ax.set_xlabel(\"capacity scenario\")\n",
    "\n",
    "ax.set_title(\"means of results\")\n",
    "ax.set_ylim([0, 2e6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of a Stochastic Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Expected costs of ignoring uncertainty (ECIU)\n",
    "\n",
    "in some literature also defined as the Value of Stochastic Solution (VSS). Can be used interchangeably.\n",
    "\n",
    "The natural question to ask is how much difference it really makes to the quality of the decisions reached if I use a stochastic problem instead of a deterministic problem?\n",
    "\n",
    "The ECIU measures the value of using a stochastic model (or the expected costs of ignoring uncertainty when using a deterministic model).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "portfolios = pd.DataFrame()\n",
    "costs = pd.Series()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the naive problem (usually -- the expected value problem (EV))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# can be anything (e.g., the 'med' scenario). A texbook way is to take expected value of uncertain parameter.\n",
    "\n",
    "naive_scenario = sum(pd.Series(gas_prices) * pd.Series(probability))\n",
    "naive_scenario\n",
    "# naive_scenario = gas_prices[\"med\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### solve naive problem (deterministic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario = \"naive\"  # naive problem (in literature often EVP for Expected Value Problem, if the naive assumption is the expected value)\n",
    "gas_price = naive_scenario\n",
    "\n",
    "n = prepare_network(cts, gas_price)\n",
    "\n",
    "n.optimize(solver_name=solver_name)\n",
    "\n",
    "portfolios[scenario] = n.generators.p_nom_opt\n",
    "costs[scenario] = n.objective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pd.set_option(\"display.precision\", 10)\n",
    "portfolios\n",
    "# costs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### solve stochastic problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario = \"SP\"  # SP for Stochastic Problem\n",
    "gas_price = gas_prices[base_scenario]\n",
    "\n",
    "n = prepare_network(cts, gas_price)\n",
    "prepare_stochastic_model(n)\n",
    "\n",
    "n.optimize.solve_model(solver_name=solver_name)\n",
    "\n",
    "portfolios[scenario] = n.generators.p_nom_opt\n",
    "costs[scenario] = n.objective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "portfolios"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve stochastic problem constrained by the naive solution "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario = \"SP-constrained\"\n",
    "\n",
    "gas_price = gas_prices[base_scenario]\n",
    "n = prepare_network(cts, gas_price)\n",
    "prepare_stochastic_model(n)\n",
    "\n",
    "n.generators.p_nom_extendable = False\n",
    "n.generators.p_nom = portfolios.loc[n.generators.index, \"naive\"]\n",
    "# n.generators.T\n",
    "\n",
    "n.optimize.solve_model(solver_name=solver_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# don't forget to add the capital costs of the (fixed) generators portfolio\n",
    "c = \"Generator\"\n",
    "ext_i = portfolios[\"naive\"].index\n",
    "cost = n.df(c)[\"capital_cost\"][ext_i]\n",
    "cost_of_portfolio = (n.generators.p_nom * cost).sum()\n",
    "n.objective += cost_of_portfolio\n",
    "n.objective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "portfolios[scenario] = (\n",
    "    n.generators.p_nom\n",
    ")  # just a fixed copy of naive problem's solution\n",
    "costs[scenario] = (\n",
    "    n.objective\n",
    ")  # must be >= than the stochastic solution's costs, because you do dispatch with the suboptimal first-stage decisions\n",
    "\n",
    "costs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute ECIU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ECIU (or VSS) in M euro\n",
    "eciu = (costs[\"SP-constrained\"] - costs[\"SP\"]) / 1e6\n",
    "# ECIU in % of stochastic solution\n",
    "eciu_pp = eciu / (costs[\"SP\"] / 1e6) * 100\n",
    "\n",
    "print(\n",
    "    f\"ECIU: {round(eciu, 3)} Meuro \\nwhich is {round(eciu_pp)}% of stochastic solution's costs\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Expected Value of Perfect Information (EVPI)\n",
    "\n",
    "If system planner knew at the first stage which scenario will play out, it could optimize an expansion plan (i.e. that results in lower cost) for that scenario.\n",
    "\n",
    "The expected value (and the corresponding mathematical problem) of such solution is denoted in the literature as „wait-and-see” solution (or wait-and-see (WS) problem).\n",
    "\n",
    "The difference between the (probability-weighted) wait-and-see solutions and the here-and-now (stochastic) solution represents the added value of information about the future (i.e., the expected profit).\n",
    "\n",
    "*modelling perspective*: How much the expected costs could be reduced if system planner in the first stage knew exactly which scenario would happen?\n",
    "\n",
    "*economic perspective*: An upper bound to the amount that should be paid for improved forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "portfolios = pd.DataFrame()\n",
    "costs = pd.Series()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve Wait-and-See problems\n",
    "where Wait-and-See (WS) is a standard textbook name for individual determinic problem (i.e. running a single scenario)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for scenario in scenarios:\n",
    "    gas_price = gas_prices[scenario]\n",
    "    n = prepare_network(cts, gas_price)\n",
    "\n",
    "    n.optimize(solver_name=solver_name)\n",
    "\n",
    "    if results is None:\n",
    "        results = pd.DataFrame(columns=n.generators.index)\n",
    "        results.index.name = \"scenario\"\n",
    "\n",
    "    portfolios[scenario] = n.generators.p_nom_opt\n",
    "    costs[scenario] = n.objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### compute the expected value of wait-and-see scenario costs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = sum(costs * pd.Series(probability))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### solve stochastic problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario = \"SP\"  # SP for Stochastic Problem\n",
    "gas_price = gas_prices[base_scenario]\n",
    "\n",
    "n = prepare_network(cts, gas_price)\n",
    "prepare_stochastic_model(n)\n",
    "\n",
    "n.optimize.solve_model(solver_name=solver_name)\n",
    "\n",
    "portfolios[scenario] = n.generators.p_nom_opt\n",
    "costs[scenario] = n.objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute EVPI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# EVPI in M euro\n",
    "evpi = (\n",
    "    costs[\"SP\"] - ws\n",
    ") / 1e6  # must be >=0 because improved information cannot make the decision maker worse\n",
    "# ECIU in % of stochastic solution\n",
    "evpi_pp = evpi / (costs[\"SP\"] / 1e6) * 100\n",
    "\n",
    "print(\n",
    "    f\"EVPI: {round(evpi, 3)} Meuro \\nwhich is {round(evpi_pp)}% of stochastic solution's costs\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing the ECIU and EVPI metrics\n",
    "\n",
    "ECIU: an investment decision is made when uncertainty is **ignored**. \n",
    "The ECIU is **the additional expected cost of assuming that future is certain**.\n",
    "\n",
    "EVPI: an investment decision is made after uncertainty is **removed**.\n",
    "The EVPI is the **expected cost of being uncertain about the future**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "",
   "language": "python",
   "name": ""
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
